Parsing mPFC: A Joint Meta-Analytic and Graph Theoretic Approach

Abstract

Valuation effects are consistently observed in medial prefrontal and posterior cingulate cortex (mPFC and PCC). The spatial extent of these effects is mostly indistinguishable from the default mode network (DMN) in existing meta-analyses. However, little is known about how valuation effects fit within the broader functional architecture of mPFC and PCC, or whether that architecture is consistent or idiosyncratic across individuals. Here we complement a meta-analysis with fMRI-based graph theoretic approaches to subdivide mPFC and PCC at the single-subject level. Our results suggest the functional topography of mPFC has substantial variability across individuals. This highlights the potential usefulness of estimating brain effects at the individual level in this region, and points to limitations of aggregative methods such as coordinate-based meta-analysis in determining whether valuation and DMN effects emerge from common or separable brain systems. Our approach shows promise in addressing this issue through future manipulations of valuation.

Introduction

Studies of decision making have consistently associated activity in ventromedial prefrontal cortex (vmPFC) with the subjective value (SV) assigned to decision outcomes (Bartra, McGuire, and Kable 2013; Clithero and Rangel 2014; Hiser and Koenigs 2018; Kable and Glimcher 2007; Levy et al. 2011). However, the distributed cortical valuation system also includes posterior cingulate (PCC) and more extended mPFC regions, which collectively show significant spatial overlap with the default mode network (DMN) along the medial wall (M. D. Fox et al. 2005; Laird et al. 2009; B. T. Yeo et al. 2011). This shared neural coverage has prompted the question of whether (and to what extent) these two systems can be dissociated.

Previous meta-analytic work has provided important insights on the psychological phenomena attributed to mPFC subregions. Independent coordinate-based meta-analyses have found mirroring spatial coverage in mPFC for valuation (Bartra, McGuire, and Kable 2013) and DMN (Laird et al. 2009). Whole brain analyses have shown that while more ventral mPFC regions are increasingly associated with reward-related elements of decision making, vmPFC strongly coactivates with DMN regions across the literature, with no clear topographical distinction between these phenomena (Vega et al. 2016). More focused meta-analytic work aimed at disentangling DMN from subjective value has determined that these systems are inseparable within mPFC (Acikalin, Gorgolewski, and Poldrack 2017). These findings highlight the multifaceted embedding of psychological constructs in subregions of mPFC, and suggest that DMN and valuation might indeed be subserved by the same system. However, recent activation-based meta-analytic work has demonstrated considerable domain specificity in this area (Kragel et al. 2018). While informative, a key limitation in distinguishing constructs through this approach is that meta-analyses often rely on information derived from group-averaged data. Averaging has traditionally been beneficial in identifying tendencies in brain function and organization when faced with short scanning sessions per subject. However, it has been shown that mPFC is subject to considerable idiosyncratic cortical folding (Zilles, Palomero-Gallagher, and Amunts 2013) and inter-subject functional variability (Mueller et al. 2013) compared to the rest of the brain, thus making group-averaging studies difficult to interpret.

Seeking to remedy these shortcomings, recent work has prescribed relevance to the analysis of single subjects in fMRI. Individual analyses of fMRI data have identified idiosyncratic, reliable, and valid functional organization that would otherwise be blurred in aggregative estimations (T. O. Laumann et al. 2015; Gordon et al. 2017; Gratton et al. 2018). Furthermore, subject-specific network arrangements have been found to predict behavioral characteristics (Kong et al. 2018). In regards to DMN, the trend of individualized analyses has led to finer idiosyncratic subdivisions of the DMN through careful selection of seed-based functional connectivity (Braga and Buckner 2017). It is thus possible that the indisitinguishable overlap of DMN and valuation effects can be attributed to a lack of spatial resolution that is better understood at the individual level. With this in mind, an important first step in disentangling these phenomena would be to determine the degree of topological heterogeneity of the DMN present within mPFC.

A persistent issue in analyzing mPFC in individuals is the various factors inducing signal dropout in its more ventral segments (Logothetis 2008). As such, there is value in examining all possible sources of covariation that compose DMN, instead of estimating topologies one seed at a time. Connectome-based analyses of resting state functional connectivity (rsFC) have been fruitful in characterizing individualized functional topologies that match task-induced activity (Gordon et al. 2017; T. O. Laumann et al. 2015; Tobyne et al. 2017). A popular approach to analyze these connectomes is to rely on graph theoretic methods, which provide an organic avenue to understand brain network dynamics (D. S. Bassett, Zurn, and Gold 2018). From these, community detection algorithms have been successfully utilized to section brain networks into cohesive substructures (Garcia et al. 2017). Such communities represent clusters of network nodes that are more connected with each other than with the rest of the network (Fortunato and Hric 2016). Among the algorithms used to estimate these communities, modularity has been widely effective in subsectioning brain networks into multiple groups (Garcia et al. 2017). However, in attempting to refine the DMN topology at the individual level, dividing these regions into what belongs and does not belong to DMN makes dissection through spectral partitioning (SP) a more valuable option. In short, spectral partitioning performs eigenvalue decompositon on an adjusted form of a correlation matrix (Laplacian matrix). The resulting eigenvector corresponding to the second-to-last eigenvalue (i.e. the Fiedler vector) provides values to divide the network into a positive and a negative community (Fiedler 1975; Belkin and Niyogi 2003). SP has recently been used to characterize the posterior-anterior functional gradient of the insula using rsfMRI data (Tian and Zalesky 2018). This method can thus capture a refined group of brain regions that more cohesively share activation patterns during rest.

In this study, we aim to subsection regions commonly attributed to both DMN and subjective value into subject-specific DMN and non-DMN partitions, quantifying the degree of topographical heterogeneity present in this network within and across individuals. This description will allow researchers to generate more precise topographic targets for future studies of decision making. We do this by capitalizing on the respective strengths of meta-analytic and subject-specific analyses of brain data. First, we define regions of interest (ROI) by identifying overlapping activation throughout the literature projected to an established brain atlas (M. F. Glasser et al. 2016). We then produce a rsFC networks of all the surface vertices within these ROIs for each individual resting state dataset from the Human Connectome Project (HCP; citation), and estimate the extent of DMN coverage through spectral partitioning. The resulting patterns are then evaluated for consistency across methods, throughout each individual’s time series, and across subjects. Finally, we show some interesting properties of SP related to the stability of brain networks.

Methods

Meta-analysis

We used data from metan-alyses that gathered peak brain coordinates of activity from studies of valuation (Bartra, McGuire, and Kable 2013) (n = 27 studies) and DMN (Laird et al. 2009) (n = 77). These represent the surviving areas post-statistical thresholding from each study. For each peak volumetric coordinate, 10 mm sphere masks were generated and then projected to a standard cortical mesh (fs_average, 32,000 vertices) using FreeSurfer. The sphere sections that overlapped with the cortex (should I add the details of the overlap?) were mapped to an atlas of the human cortical surface (360 regions)(M. F. Glasser et al. 2016). This produced a list of standardized parcels that were reported on each study.

With these lists, we identified the areas that were more consistently shared between DMN and valuation. For each literature, we summed the number of times each parcel was reported across studies, and performed a one-way chi-squared test of proportions on each parcel with a null probability equivalent to the bootstrapped (5000 iterations) mean of each literature’s report counts. In other words, we checked which regions were significantly (p < 0.05) more reported in their respective literatures than their counterparts. Regions of interest were then determined based on which of the remaining regions were common across literatures.

Next, we identified areas that were solely related to either literature. We limited the search space to the parcels that were contained within the DMN parcellation defined by Yeo et al. (2011). We permuted the label of every study (i.e. DMN or valuation) 5000 times while preseving the total number of studies observed in each domain, and on each iteration stored the maximum statistic from a ROI-wise chi-squared test of proportions on the number of prmuted studies that had reported each ROI. This gave us a null distribution of chi-squared values that was used to determine which ROIs were more significantly prevalent on each literature.

fMRI Data

In order to quantify the intrinsic connectivity of these ROIs in the brain, we acquired resting-state fMRI data from the Human Connectome Project (HCP) (Van Essen et al. 2012) Q6 release (n = 80, randomly sampled from the total pool of 469 available subjects). Each subject’s data was acquired over two days at Washington University in St. Louis on a Siemens CONNECTOM Skyra MRI scanner (Siemens, Erlangen, Germany). Four resting state sessions (TR = 0.720 s, TE = 33.1 ms, FA = 52°, multiband factor = 8, 72 slices, 2 mm isotropic voxels) were each comprised of 1200 TRs for a total scan time of 14 min 33 s, with left-right and left-right phase encoding implemented on each day. We used this scanning structure to estimate each subject’s topology based on session, day, and overall data (4800 TRs).

Beyond the in-house minimally preprocessed pipeline from HCP (M. F. Glasser et al. 2013), which includes considerable motion correction, temporal denoising, highpass filtering (0.0005 Hz threshold), and MNI152-based normalization, all scans went through a number of additional refinements (Tobyne et al. 2017). These included band-pass filtering (allowed frequencies ranged from 0.009 and 0.08 Hz), as well as mean greyordinate signal regression (Burgess et al. 2016). Only subjects with both left-right and right-left phase encoding were included (i.e. subjects with four rsfMRI sessions). In addition, only datasets with either low motion levels (under 1.5mm) or under 0.5 mm mean framewise displacement (FD) were used. Volumes that displayed an FD of over 0.5 mm were considered as spikes and removed from the subject’s dataset, and subjects whose spike counts were above 5% of their total data were excluded from analyses. Finally, data acquired on the same day (i.e. left-right and right-left phase encoding session) were temporally demeaned. Each subject’s brain was comprised of 32k vertices per hemisphere. We retained only the cortical surfaces, which resulted in 59,412 total surface vertices per subject.

Within-ROI Network Definition

To establish each subject’s network, we selected all the vertices contained within the ROIs (n = 5,319 per subject) and computed the Pearson correlation of the time series for every pair of vertices. This produced a weighted network for each subject, where surface vertices were the nodes, and edges the correlations among all of them. Next, we applied Fisher’s r to z transformation, and took the exponential of the correlations so that all weights were positive while maintaining the shape of the original correlation distribution, and retained the edges whose adjusted p-values remained significant after correction for multiple tests (FDR < 0.05). As part of the evaluation step, this procedure (and the community detection outlined below) was applied at four levels: 1) each session separately (1200 TRs each); 2) the concatenated timeseries from each pair of daily sessions (2400 TRs total); 3) each subject’s full dataset (4800 TRs total); and 4) on each step of a sliding window analysis (see Evaluation for more details).

Community Detection

Communities (i.e. clusters) were identified by means of spectral partitioning (Fiedler 1975; Belkin and Niyogi 2003). First, each adjusted network was represented as an \(n\) x \(n\) correlation matrix (where \(n\) = number of vertices within the ROIs, 5,319). This matrix was then transformed into its normalized Laplacian form

\[\begin{aligned} L = I - D^{-\frac{1}{2}}WD^{-\frac{1}{2}} \end{aligned}\]

Where \(I\) is the identity matrix, \(D\) is a diagonal matrix containing the strength of each vertex (i.e. the sum of its correlations with every other vertex), and \(W\) is the correlation matrix. We then computed the eigenvalues and eigenvectors of the normalized Laplacian, and used the eigenvector associated with the second-to-lowest eigenvalue to divide the network into two. This vector (from now on called the ‘Fiedler Vector’) provides a set of positive and negative values to binarize the network with, and guarantees that the resulting communities are connected (Fiedler 1975). Given that this data-driven method is agnostic to the prevalent function of each community (i.e. DMN or not DMN), we defined the DMN community as the one containing the one DMN-specific area from the meta-analysis (ensuring that the labeling of the DMN community was consistent across subjects). Importantly, given the high density of these networks (see results), spectral partitioning was unlikely to face the issues associated with its use in sparse networks (Fortunato and Hric 2016).

In order to evaluate the validity of the resulting partitions, we also estimated these communities using the more traditional approach of modularity maximization using the algorithm from Clauset et al. (2004). This method maximizes a quality function that compares the strength of the connection between any two vertices against a null model of their probability of being randomly connected to any other vertex in the network. The method heuristically iterates through many possible combination of vertices, and selects the clustering that maximizes the network’s modularity. Thus, unlike SP, modularity can fractionate a network into many functional groups. If the partitions from the bounded and unbounded methods show high level of agreement, we can gain confidence in that we are indeed fractioning the ROIs into those belonging to DMN and not.

Partition Evaluation

In order to evaluate the stability and topographical heterogeneity of the partitions within and across subjects, we chose the adjusted rand index (ARI) (Hubert and Arabie 1985) as the reference metric of agreement between two clusterings. The adjusted rand index estimates the rate at which every pair of nodes in partition A were equally assigned to the same/different parcels on partition B, over all possible changes (e.g. clustered together on A, but separately on B). The resulting ratio is compared to a baseline given by the expectation based on a random vertex assignment for an equal number of clusters across partitions. Because of this, the adjusted rand index is a simple metric that determines the percentage of agreement between any two partitions that is agnostic to the labeling scheme. (SHOULD WE INCLUDE THIS? The mutual information results would probably need to go in the supplemental materials) On the other hand, variation of information gives a general sense of the amount of information to be gained by the complementary community partitioning (Meila, 2007), measured by the addition of the conditional entropies per combination of clusters across partitions. While the rand index is more easily interpreted, it suffers from the assumption that the baseline should have an equal number of vertices per cluster across partitions, which is not necessarily the case. Therefore, I used variation of information as a more robust similarity metric (Fortunato & Hric, 2016).

We performed a number of comparisons among partitions. First, we computed the degree of agreement between SP and modularity per subject to validate the resulting functional patterns. Modularity and SP show overfitting and underfitting tendencies, respectively, in their community detection performance in a diverse number of network types (Ghasemian, Hosseinmardi, and Clauset 2018). Their alignment would thus give us confidence in the validity of the resulting topologies. Next, we compared the functional organization derived from the whole time series among all subjects, and report the mean RI for the group. We then estimate the level of agreement between scans taken across days (i.e. pairs of daily left-right and right-left phase encoding sessions). If the functional organization estimated by SP is indeed subject-specific, we should see higher agreement for the same subject across days than either day compared to every other subject. We then extend this idea by computing the overall agreement across all sessions. Similarly, we expect that sessions within subject will show higher agreement than between subjects. Finally, we perform a sliding window analysis (20 min windows, 1 min increments, mean number of windows per subject = 37) 8) comparing each window’s functional pattern against the clustering derived from the whole subject’s data. This approach should give us a sense of the dynamics of the functional organization of mPFC during rest. We complement this analysis by calculating the proportion of times each surface vertex was affiliated with the DMN community across time, as well as computing the flexibility of each node (i.e. the number of times each vertex changed affiliation across time). Given that the functional heterogeneity of brain areas is unevenly distributed across the brain (Mueller et al. 2013), we also performed these comparisons for PCC and mPFC separately. Given what we know about these general regions, we expected PCC to be well aligned across participants, and stable within subjects. On the other hand, we hypothesized low agreement of mPFC partitions across participants, but with considerable stability across days and sessions within subject. We compared the agreement within subjects against across subjects by computing the ratio of the mean RI within and between subjects. Ratios close to one denote good alignment across subjects, while ratios significantly higher than one suggest that partitions per subject were more stable than what you would expect against every other subject.

Given the descriptive relevance of Fiedler eigenvector values in previous fMRI research (Tian and Zalesky 2018), we tested the association of their absolute magnitude with the proportion of times a vertex was part of the DMN community during the sliding window analysis. To do this, we performed a mixed effects logistic regression model in which each subject’s eigenmap values were used as random predictors of the proportion of times they participated in the DMN community.

Estimation of Noise Effects

BOLD signal recording through fMRI in mPFC in general, and vmPFC in particular, suffers from technical issues that promote signal dropout (Logothetis 2008). It was thus important to quantify the potential influence of noise in our findings. We examined this through two complementary steps. First, we computed the temporal signal-to-noise ratio (tSNR) of all the vertices in the ROIs from HCP minimally preprocessed data (i.e. prior to the complementary denoising described in the “fMRI Data” section). For each subject, we compared the mean tSNR between mPFC and PCC through a permutation analysis (5000 iterations), and report the percentage of subjects who had a significantly lower mean tSNR in mPFC. Next, we estimated which vertices within mPFC had a stronger relationship to the overall dynamics of DMN. To do this, we computed the mean time series of each subject’s DMN map, which was then predicted by the time series each mPFC vertex in a mass univariate model. The resulting \(R^2\) values from each model were extracted to determine the total variance of the mean DMN time series explained by each individual vertex. If high noise in more ventral mPFC areas is affecting our results, we should see a gradient of \(R^2\) values that decreases ventrally across subjects. We then compare this dorsal-to-ventral gradient to one displaying tSNR values, as the latter should convey a linear relationship between noise and dorsal-to-ventral topographical placement.

Seed-based FC vs Community Detection

A potential downfall of the present approach is that it can be computationally demanding, while more direct seed-based functional connectivity methods could provide sufficient approximations of an individual’s DMN orgnization. We examined this issue by estimating mPFC functional maps from its vertex-wise functional coupling with the mean PCC activation. For each subject’s concatenated day 1 and day 2 scans, we computed the mean time series of all vertices contained within the PCC region of a standard DMN arrangement (B. T. Yeo et al. 2011), and correlated it (Pearons) with the time series of every mPFC vertex from our ROIs. The resulting day 1 and day 2 mPFC maps were then correlated within and between subjects in three ways: 1) between correlation-based maps; 2) between community-based maps; and 3) between correlation- and community-based maps. These spatial correlations are aimed to determine the day-to-day stability of each approach, as well as their level of agreement. These methods should produce somewhat similar arrangements, but the one displaying the highest within-subject agreement across days should be preferred (for example, see (Kong et al. 2018; Mueller et al. 2013)). We thus compared the mean within-subject inter-day correlation coefficients of each method through a paired permutation analysis (5000 iterations).

Anatomical Correspondence

Finally, we examined whether each idiosyncrasies in functional topology corresponded to outstanding anatomical features. To do this, we correlated each subject’s Fiedler vector values with their respective cortical thickness, cortical folding, and myelin density (quantified as the ratio of T1 to T2 magnitudes at every given vertex). These simple correlations are aimed to provide a coarse sense of this possible relationship, which can be studied in depth in the future.

Results

Meta-analysis

We performed a coordinate-based meta-analysis to limit the network search space to a set of overlapping and unique standard regions across DMN and SV literatures (see Methods). Volumetric coordinates from a set of DMN and SV studies were projected onto a cortical surface, and mapped to a multimodal cortical parcellation (M. F. Glasser et al. 2016) to produce a list of brain areas reported per study. Using each literature’s mean number of reports per area as the null probability (DMN = 12, 95% CI [11.53, 12.92]; SV = 5, 95% CI [4.56, 5.18]), one-way chi-squared tests of proportions were performed to identify regions that were significantly present in their respective literatures. The resulting lists were used to identify the set of common brain areas. Unique areas were then determined by permuting the affiliation of each study (i.e. DMN or SV) to create a null distribution of chi-squared statistics, which as used to identify regions more often reported in either literature. Importantly, we chose to examine these areas bilaterally, even if their significance was present in only one hemisphere.

Figure 1 shows the resulting ROIs from these analyses. As expected, overlapping areas included 15 bilateral regions in PCC (31a, 31pv, RSC, d23ab) and mPFC (9m, 10r, 10v, OFC, a24, d32, p24, p32, s32), as well as a region in IPL (PF) and TPJ (PGi) (all p < 0.05). A region in visual cortex (V2) was also significant, but was removed due to its know unrelated function. Area 7m was the only one to show unique affiliation, which was significant for DMN (\(\chi^2\) = 9.12, p < 0.001).

Figure 1. Template atlas from Glasser et al. (2016) (top), as well as the common and unique brain regions associated with DMN and subjective value (bottom). Overlapping ROIs (in blue) were defined as the common brain regions that were significantly present in their respective literatures, while a single DMN-specific region (7m, in red) was identified from a permutation analysis comparing these literatures.

Partition Comparison: Full Datasets

With these 16 ROIs (15 overlapping, 1 DMN-specific), we estimated each individual’s DMN coverage using the whole time series (4800 TRs). Each subject’s network was constructed by performing pairwise correlations among the time series of all the vertices contained within the ROIs. After correcting for multiple tests, we performed spectral partitioning and modularity on these networks to identify cohesive functional communities denoting DMN and non-DMN areas. On average, networks had high densities (i.e. highly connected) regardless of the amount of data to produce them (mean = 0.96, SD = 0.0099736). This degree of density is important, as these methods are known to face challenges when computed on sparse networks (Fortunato and Hric 2016). Figure 2 shows a stereotypical subdivision of the ROIs for a single subject (100307) using spectral partitioning.

Figure 2. Brain partition for an example subject (100307). Fiedler vector values (top) are mapped onto the brain surface, dividing it into positive and negative communities. The bottom brain shows the binarized Fiedler vector, with yellow areas denoting the DMN community (as indicated by coverage of area 7m).

This first step just gives an idea of the level of agreement within and between subjects. As a first measure, we use the adjusted rand index (RI) to estimate the level of agreement between spectral partitioning and modularity. RI measures the proportion of pairs of nodes in a network that were either clustered together or separately in both partitions, over all possible affiliation changes (e.g. clustered together in partition 1, but separately in 2). It is worth noting that the adjusted rand index has well known limitations (Fortunato & Hric, 2016). To compensate for this, we also performed all the upcoming analyses with the variation of information index (VI). This distance measure denotes how much information is gained by adding the second partition (thus, lower values are better), and is known for its robustness. Both RI and VI analyses agreed in the results they yielded, but we only report RI due to its easier interpretability. This first plot shows the distribution of RI values between modularity and spectral partitioning for each subject.

The total number of subjects is 100. The distribution suggests that both clustering methods were in high agreement across subjects (mean = 0.85 sd = 0.11). The importance of this results lies in that while spectral partitioning forces a bipartition of the network, modularity does not. This adds evidence to the division of these areas into DMN and non-DMN during rest.

Next, we examine how similar the partitions for PCC and mPFC are across subjects. Qualitative inspection of the brain community topologies shows good alignment for PCC, with a topologically similar mPFC pattern that shifts topographically across subjects. The following similarity matrix shows RIs for every pairwise comparison among all subjects for PCC (lower matrix) and PFC (upper) separately. Each row/column element is a single subject.

The results match the qualitative observations mentioned above, with overall higher agreement for PCC (paired permutation, p = 0; Cohen’s D = 4.21). However, the heterogeneity of mPFC could arise from noise in the signal. As a way to test for subject specificity, we used each subject’s session data to compare day 1 vs day 2 clusters, as well as a session-based analysis.

Partition Comparison: Day 1 vs Day 2

HCP resting state data is separated into 4 runs (2 per day of scanning). This section explores how much agreement there is between partitions computed on each day. The similarity matrix on the left shows the pairwise level of agreement for the overall data among the example subjects, while the right one divides this data into PCC and mPFC.

So far both visualizations show considerable subject specificity, and the regional division agrees with the whole-data similarity matrix. Importantly, this indicates that both PCC and mPFC partitions were more similar within- than between-subjects (e.g. RI values higher for each 2x2 diagonal submatrix than the off-diagonals). Now, how do we check whether the within subject consistency is significantly higher than across subjects? One way is to divide each subject’s RI-agreement between days by the RI values with the days from the other subjects. This gives an RI ratio per subject that denotes the relative specificity of their partition. That means, partitions that align well across subjects should yield ratios close to 1, while subject-consistent and unique partitions would produce high ratios.

As expected, we see that while both regions had good specificity (ratios > 1, denoted by the horizontal line), the ratios for PCC were significantly smaller than mPFC for all subjects (permutation, p = 0; Cohen’s D = 2.9343989). However, the denominator is given by a single RI. If these topological patterns are indeed subject-specific for mPFC, these results might generalize to session specific data. This would give us a better estimate of within-subject RI variance to place in the denominator.

Partition Comparison: Session

Next, each session is evaluated independently. The similarity matrix on the left shows how each session’s data corresponded within and between subjects for all ROIs combined, while the one on the right shows PCC and mPFC separately. It is worth noting that each session is ~14 mins long, which is below the stability threshold for fMRI based modularity (Gordon et al., 2017). However, we also perform a sliding window analysis below that encompasses 20 mins of data per window.

As before, we see higher agreement levels within subjects (with some exceptions). Simple permutation comparisons confirmed this intuition for the all-ROI (p = 0; Cohen’s D 2.01), PCC (p = 0; Cohen’s D = 1.59), and mPFC (p = 0; Cohen’s D = 1.85) divisions. In order to quantify the consistency and specificity of these partitions, we again applied the RI ratio index, but this time using the mean RI among all of a subject’s sessions in the denominator. The following plot shows the results for PCC and mPFC separately.

We can see again that the ratios for mPFC were higher than those of PCC (permutation, p = 0; Cohen’s D = 2.66), but both were above 1. This suggests that the PCC is well aligned across subjects (as mentioned before), while mPFC is consistent within subject, but quite different from other subjects.

Another way to determine the heterogeneity of these results is based on Kragel et al. (2018), who quantified the dissimilarity of brain patterns for behavioral domain, subdomain, and study (or how generalizeable the patterns are at each hierarchical level). In our case, for 100 subjects 300 matrices are created (200 for daily comparisons, 100 for within subject comparisons). Each matrix is binarized to compare a subject’s level of similarity within their daily data vs across every other subject’s sessions. The 300 matrices are vectorized and used as predictors in a GLM with a quasibinomial distribution, where the dependent data is the full vectorized RI similarity matrix. This is a simplified version of a fully-dimensional model that considers all possible comparisons at the session level, which would otherwise have S + (S x (S-1)/2) parameters (with S = session). For this model, a significantly negative subject predictor implies that the odds of the partitions being similar across subjects are much lower than those within each subject (across sessions). Thus, subject specificity would dictate that all subject-level coefficients should be significant and negative.

The plot above shows the estimated coefficients from the linear model for day and subjects (R-squared for the model = 0.53). As we would expect given the subject specificity of these partitions, 0.99 percent of the subject coefficients were significant (a one sample t-test confirmed this notion). Conversely, 0.38 of the coefficients for daily comparisons were significant (non-significant t-test). These results suggest that the partitions were better explained as being subject-specific. That said, the current version of the model is fit with the data from all ROIs, and running it separately for PCC and mPFC might be fruitful in the future.

In order to diagnose the model, we checked the independency of observations through the rank of the design matrix, which was 300 (full rank for 100 subjects should be 300) for the present analysis. Further diagnostics included a qqnorm check on the model residuals (result: mostly normal), and variance inflation factors (VIF, also used by Kragel et al., 2018). The VIF is a measure of colinearity of the regressors, and it “shows the degree to which variance in each parameter estimate is increased due to partial colinearity with linear combinations of other regressors”. While low (<2) VIF values are preferable, shared variance is unavoidable given the hierarchical structure of the model. Regardless, Kragel et al. take it as a good sign that their VIF are finite (?), which is also the case here (mean = 40.6185275). That said, these are above 10, which is a standard upper threshold (potential downside). We will have to reconsider this portion.

Stability per subject

It is also worth considering how stable these measures are, so we performed a sliding window analysis per subject (20 min windows shifting at ~1 min throughout the whole ~ 1hr time series). The first two plots show how well each window’s clustering agreed with the partition derived from the whole data (in RIs, left), as well as how variable (in standard deviations) this index was throughout each time series (right).

While some subjects show considerable variability, the overall trend is of fair stability. An interesting feature of spectral partitioning is that the Fiedler vector, which provides signed values to divide a network into positive and negative communities, was tightly related to the proportion of times any given cortical vertex was affiliated with the DMN throughout the fMRI sessions. Meaning, the higher the absolute value in the Fiedler vector for a given vertex, the more persistent its relationship was with its corresponding community in time. The following plot shows this relationship for each subject (gray smoothed lines), as well as the fitted values from a logistic mixed effects model (black line; model: prop. DMN ~ Fiedler vector, subjects as random intercepts; beta = 757.38, SE = 1.68, p < 0.001).

With the exception of some noisy subjects, the sigmoidal functions show a tight relationship between the two. This brings attention to the possibility of using the Fiedler vector as a surrogate for stability, thus reducing the computational needs for future analyses.

While these results show promise, so far we have only looked at stability for all the ROIs. The next three plots show similar data, but for PCC and mPFC separately. The middle plot compares the subject-wise mean RI for PCC and mPFC. While mPFC shows an overall significantly lower level of agreement with the overall partition, the variance throughout the time series is small.

This suggests that, while the agreement of each window with the overall time series was stable for both PCC and mPFC (low SD, right plot), the overall agreement of mPFC is lower (paired permutation, p = 0; Cohen’s D = 3.67). This result could be due to signal degeneration in mPFC (that said, it is worth noting that the topological patterns remain mostly the same for both PCC and mPFC throughout the SW analysis). To account for this, we computed the temporal signal to noise ratio (tSNR) for each vertex in each subject’s minimally preprocessed fMRI data (i.e. mean of each vertex’s time points over their standard deviation). The next plot shows the pooled subject-wise mean tSNR for PCC and mPFC vertices, and indicates that the tSNR is lower for mPFC than PCC, which could negatively affect our results.

To formally examine this, for each subject we compared the tSNR values of all vertices between PCC and mPFC (100 permutation tests, 5000 iterations). About 99 percent of the comparisons proved significant at alpha = 0.05 (FDR corrected; mean Cohen’s D = 0.83). In addition, we found that the degree of signal to noise for each vertex showed significant negative correlations with flexibility, a graph-theoretic index that measures the proportion of times a node in a network changed affiliation across partitions (mean correlation across subjects = -0.23). This suggests that, for any given subject, vertices with low SNR resulted in variable participation between the communities, and could explain why mPFC had lower stability levels.

To further quantify the potential effect of noise in our results, we computed the variance explained by the mean time series of the DMN community on each DMN vertex separately through a simple linear model. Interestingly, there is a strong linear relationship between the Fiedler vector (the divisive vector from spectral partitioning) and the R-squared values (these match perfectly on the brain when inspected qualitatively). The following plot shows this relationship, with each line denoting a subject.

Next, we compared the R-squared values for PCC and mPFC separately per subject. The following plot shows that for each subject, R-squared values tend to be higher for PCC (100 percent of Wilcoxon ranksum tests across subjects had p values < 0.001). This is also true for Fiedler vector values, but to a lesser extent. These differences had moderately high effect sizes (mean Cohen’s D = 0.98). These results suggest that overall, vertices belonging to the mPFC DMN community are less in synchrony with the activity of the community as a whole than those in PCC.

This begs the question of whether there is a “gradient of noise” along the z-axis on mPFC (vmPFC explaining less variance?). The plots below show the density of R-squared (left) and SNR (right) values for every mPFC vertices across all participants. The R-squared values have no clear relationship with the ventral-dorsal axis (with some low values grouped in vmPFC). However, the coefficient of variation derived from minimally preprocessed data does show a trend (namely, the more ventral, the lower the SNR). This is notable because the R-squared values are estimated from data with added noise reduction steps. These results contextualize the degree of synchrony of mPFC vertices in terms of potential noise. While still an important factor, it is unlikely that our stability results could be due mostly to highly structured mPFC noise.

Correlation vs community detection in mPFC

Next, it is worth examining what we gain from using graph theoretic approaches over a simple seed-based functional connectivity analysis from PCC. For this purpose, we computed the mean time series of all the vertices contained in the DMN PCC ROI from Yeo et al. (2011) parcellation, and correlated it with the activity of every mPFC vertex from our ROIs. The similarity matrix below shows Day 1 vs Day 2 comparisons for mPFC partitions generated with seed-based correlations (top-left quadrant) vs community detection (bottom-right) across all participants, as well as how much did these methods agree with each other (not much, but enough, which is good). The important features here are that both approaches are very consistent, but community detection seems slightly more stable across days (quantified in the next section). Also, the cross-method agreement is fair, meaning that both approaches are showing similar topological features, but that these are more stable using graph theoretic approaches.

The previous plot shows the inter-methodological correspondence with an example subjects. The final plot shows the spatial correlation between days for each method across all participants. Community reflects the data from the lower right quadrant, correlation belongs to the upper left, and both is the cross-approach correspondence (top right). Community detection had a high and significantly higher RI across days than a simple correlation (permutation p-value < 0.001; Cohen’s D = 1.14). This is important because it shows that we can derive very stable topological features using community detection, even though mPFC has relatively low signal to noise ratio. Low intra-subject and high inter-subject variability of DMN partitions is in agreement with previous work [T. O. Laumann et al. (2015); Gordon et al. (2017); Braga and Buckner (2017); Kong et al. (2018)). Graph-based community detection might thus be a better way to detect stable functional topologies, since high stability throughout days is a consistent feature of DMN topologies.

Structural correspondence

Finally, we examined whether these topological features were defined by the structure of the brain by correlating each subject’s Fiedler vector with their structural MRI data. We found that the Fiedler vector is uncorrelated with myelin density (mean correlation = -0.01; SD = 0.06), cortical folding (mean correlation = -0.02; SD = 0.04), and surface thickness (mean correlation = 0.11; SD = 0.07). While these results suggest that each subject’s functional topology is mostly independent of brain structure, we will perform more detailed analyses of this in the future (i.e. are there any outstanding structural features in a subject that match their unique topology?).

Discussion

For clinical implications, see Hiser & Koenigs 2018. In short, being able to identify functional idiosyncrasies without much data can yield important insights in diagnosis (meh, we always say that..) of a variety of psychological diseases associated with mPFC. It would be more interesting to identify the repercussions to damage in mPFC per individual. They mention that subgenual vmPFC is a common target for deep brain stimulation (or TMS) as a depression treatment. The refinement of the technique is thus coupled to being able to identify personalized functional topologies. They list a number of results related to mPFC functional connectivity too. In the end, their suggestions for clinical implications of mPFC studies rely on the individual measure of functional expanse. (note that social-relevant areas in their meta analysis look similar to our patterns, and the emotion part to Barbas latest)

It will also be good to note that there is concurrent debate on the ground truth of community detection. In other fields, identifying a community structure is not necessarily associated with the true organization of the network. What happens is that you might find a tangential grouping that gives you information different from your original goal (Clauset, Newman, and Moore 2004)

Mention the lack of interhemispheric distinction, citing Mackey and Petrides in regards to L-R diffs.

Acknowledgements

Data were provided [in part] by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.

References

Acikalin, M. Yavuz, Krzysztof J. Gorgolewski, and Russell A. Poldrack. 2017. “A coordinate-based meta-analysis of overlaps in regional specialization and functional connectivity across subjective value and default mode networks.” Frontiers in Neuroscience 11 (JAN): 1–11. doi:10.3389/fnins.2017.00001.

Bartra, Oscar, Joseph T. McGuire, and Joseph W. Kable. 2013. “The valuation system: A coordinate-based meta-analysis of BOLD fMRI experiments examining neural correlates of subjective value.” NeuroImage 76. Elsevier Inc.: 412–27. doi:10.1016/j.neuroimage.2013.02.063.

Bassett, Danielle S, Perry Zurn, and Joshua I Gold. 2018. “On the nature and use of models in network neuroscience.” Nature Reviews Neuroscience. Springer US, 1–13. doi:10.1038/s41583-018-0038-8.

Belkin, Mikhail, and Partha Niyogi. 2003. “Laplacian eigenmaps for dimensionalitu reduction and data representation” 1396: 1373–96.

Braga, Rodrigo M., and Randy L. Buckner. 2017. “Parallel Interdigitated Distributed Networks within the Individual Estimated by Intrinsic Functional Connectivity.” Neuron 95 (2). Elsevier Inc.: 457–471.e5. doi:10.1016/j.neuron.2017.06.038.

Burgess, G. C., S. Kandala, D. Nolan, T. O. Laumann, J. D. Power, B. Adeyemo, M.P. Harms, S.E. Petersen, and D. M. & Barch. 2016. “Evaluation of Denoising Strategies to Address Motion-Correlated Artifacts in Resting-State Functional Magnetic Resonance Imaging Data from the Human Connectome Project.” Brain Connectivity 6 (9): 669–80.

Clauset, Aaron, M. E. J. Newman, and Cristopher Moore. 2004. “Finding community structure in very large networks” 066111: 1–6. doi:10.1103/PhysRevE.70.066111.

Clithero, John A., and Antonio Rangel. 2014. “Informatic parcellation of the network involved in the computation of subjective value.” Social Cognitive and Affective Neuroscience 9 (9): 1289–1302. doi:10.1093/scan/nst106.

Fiedler, Miroslav. 1975. “A Property of Eigenvectors of Nonnegative Symmetric Matrices and its Application to Graph Theory.” Czechoslovak Mathematical Journal 25 (100): 619–33.

Fortunato, Santo, and Darko Hric. 2016. “Community detection in networks: A user guide.” Physics Reports 659. Elsevier B.V.: 1–44. doi:10.1016/j.physrep.2016.09.002.

Fox, M. D., A. Z. Snyder, J. L. Vincent, M. Corbetta, D. C. Van Essen, and M. E. Raichle. 2005. “From The Cover: The human brain is intrinsically organized into dynamic, anticorrelated functional networks.” Proceedings of the National Academy of Sciences 102 (27): 9673–8. doi:10.1073/pnas.0504136102.

Garcia, Javier O., Arian Ashourvan, Sarah Muldoon, Jean M. Vettel, and Danielle S. Bassett. 2017. “Applications of community detection techniques to brain graphs: Algorithmic considerations and implications for neural function.” bioRxiv, 209429. doi:10.1101/209429.

Ghasemian, Amir, Homa Hosseinmardi, and Aaron Clauset. 2018. “Evaluating Overfit and Underfit in Models of Network Community Structure,” 1–17. doi:arXiv:1802.10582v2.

Glasser, Matthew F., Timothy S. Coalson, Emma C. Robinson, Carl D. Hacker, John Harwell, Essa Yacoub, Kamil Ugurbil, et al. 2016. “A multi-modal parcellation of human cerebral cortex.” Nature 536 (7615). Nature Publishing Group: 171–78. doi:10.1038/nature18933.

Glasser, Matthew F., Stamatios N. Sotiropoulos, J. Anthony Wilson, Timothy S. Coalson, Bruce Fischl, Jesper L. Andersson, Junqian Xu, et al. 2013. “The minimal preprocessing pipelines for the Human Connectome Project.” NeuroImage 80. Elsevier Inc.: 105–24. doi:10.1016/j.neuroimage.2013.04.127.

Gordon, Evan M., Timothy O. Laumann, Adrian W. Gilmore, Dillan J. Newbold, Deanna J. Greene, Jeffrey J. Berg, Mario Ortega, et al. 2017. “Precision Functional Mapping of Individual Human Brains.” Neuron 95 (4). Elsevier Inc.: 791–807.e7. doi:10.1016/j.neuron.2017.07.011.

Gratton, Caterina, Timothy O. Laumann, Ashley N. Nielsen, Deanna J. Greene, Evan M. Gordon, Adrian W. Gilmore, Steven M. Nelson, et al. 2018. “Functional Brain Networks Are Dominated by Stable Group and Individual Factors, Not Cognitive or Daily Variation.” Neuron. Elsevier Inc., 439–52. doi:10.1016/j.neuron.2018.03.035.

Hiser, Jaryd, and Michael Koenigs. 2018. “The Multifaceted Role of the Ventromedial Prefrontal Cortex in Emotion, Decision Making, Social Cognition, and Psychopathology.” Biological Psychiatry 83 (8). Elsevier Inc: 638–47. doi:10.1016/j.biopsych.2017.10.030.

Hubert, Lawrence, and Phipps Arabie. 1985. “Comparing partitions.” Journal of Classification 2 (1): 193–218. doi:10.1007/BF01908075.

Kable, Joseph W., and Paul W. Glimcher. 2007. “The neural correlates of subjective value during intertemporal choice.” Nature Neuroscience 10 (12): 1625–33.

Kong, Ru, Jingwei Li, Csaba Orban, Mert Rory Sabuncu, Hesheng Liu, Alexander Schaefer, Nanbo Sun, et al. 2018. “Spatial Topography of Individual-Specific Cortical Networks Predicts Human Cognition, Personality and Emotion.” bioRxiv 36 (24): 213041. doi:10.1523/JNEUROSCI.4402-15.2016.

Kragel, Philip A., Michiko Kano, Lukas Van Oudenhove, Huynh Giao Ly, Patrick Dupont, Amandine Rubio, Chantal Delon-Martin, et al. 2018. “Generalizable representations of pain, cognitive control, and negative emotion in medial frontal cortex.” Nature Neuroscience. Springer US, 1. doi:10.1038/s41593-017-0051-7.

Laird, A. R., S. B. Eickhoff, K. Li, D. A. Robin, D. C. Glahn, and P. T. Fox. 2009. “Investigating the Functional Heterogeneity of the Default Mode Network Using Coordinate-Based Meta-Analytic Modeling.” Journal of Neuroscience 29 (46): 14496–14505. doi:10.1523/JNEUROSCI.4004-09.2009.

Laumann, Timothy O., Evan M. Gordon, Babatunde Adeyemo, Abraham Z. Snyder, Sung Jun Joo, Mei Yen Chen, Adrian W. Gilmore, et al. 2015. “Functional System and Areal Organization of a Highly Sampled Individual Human Brain.” Neuron 87 (3). Elsevier Inc.: 658–71. doi:10.1016/j.neuron.2015.06.037.

Levy, I., S. C. Lazzaro, R. B. Rutledge, and P. W. Glimcher. 2011. “Choice from Non-Choice: Predicting Consumer Preferences from Blood Oxygenation Level-Dependent Signals Obtained during Passive Viewing.” Journal of Neuroscience 31 (1): 118–25. doi:10.1523/JNEUROSCI.3214-10.2011.

Logothetis, Nikos K. 2008. “What we can do and what we cannot do with fMRI.” doi:10.1038/nature06976.

Mueller, Sophia, Danhong Wang, Michael D. Fox, B. T.Thomas Yeo, Jorge Sepulcre, Mert R. Sabuncu, Rebecca Shafee, Jie Lu, and Hesheng Liu. 2013. “Individual Variability in Functional Connectivity Architecture of the Human Brain.” Neuron 77 (3). Elsevier Inc.: 586–95. doi:10.1016/j.neuron.2012.12.028.

Tian, Ye, and Andrew Zalesky. 2018. “Characterizing the functional connectivity diversity of the insula cortex: Subregions, diversity curves and behavior.” NeuroImage. Elsevier Inc. doi:10.1016/j.neuroimage.2018.08.055.

Tobyne, Sean M., David E. Osher, Samantha W. Michalka, and David C. Somers. 2017. “Sensory-biased attention networks in human lateral frontal cortex revealed by intrinsic functional connectivity.” NeuroImage 162 (August). Elsevier Ltd: 362–72. doi:10.1016/j.neuroimage.2017.08.020.

Van Essen, D. C., K. Ugurbil, E. Auerbach, D. Barch, T. E.J. Behrens, R. Bucholz, A. Chang, et al. 2012. “The Human Connectome Project: A data acquisition perspective.” NeuroImage 62 (4): 2222–31. doi:10.1016/j.neuroimage.2012.02.018.

Vega, Alejandro de la, Luke J. Chang, Marie T. Banich, Tor D. Wager, and Tal Yarkoni. 2016. “Large-Scale Meta-Analysis of Human Medial Frontal Cortex Reveals Tripartite Functional Organization.” The Journal of Neuroscience 36 (24): 6553–62. doi:10.1523/JNEUROSCI.4402-15.2016.

Yeo, B.T. Thomas, Fenna M. Krienen, Jorge Sepulcre, Mert R. Sabuncu, D. Lashkari, Marisa Hollinshead, Joshua L. Roffman, et al. 2011. “The organization of the human cerebral cortex estimated by intrinsic functional connectivity.” Journal of Neurophysiology 106: 1125–65. doi:10.1152/jn.00338.2011.

Zilles, Karl, Nicola Palomero-Gallagher, and Katrin Amunts. 2013. “Development of cortical folding during evolution and ontogeny.” Trends in Neurosciences 36 (5). Elsevier Ltd: 275–84. doi:10.1016/j.tins.2013.01.006.